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Abstract 

We consider domination analysis of approximation algorithms for the bipartite boolean 
quadratic programming problem (BBQP) with m + n variables. A closed form formula is de- 
, veloped to compute the average objective function value s? of all solutions in 0(mn) time. 

However, computing the median objective function value of the solutions is shown to be NP- 
hard. Also, we show that any solution with objective function value no worse than si dominates 
at least 2 m+,l ~ 2 solutions and this bound is the best possible. Further, we show that such a 
solution can be identified in 0(mn) time and hence the dominance ratio of this algorithm is at 
least \. We then show that for any fixed rational number a > 1, no polynomial time approxima- 
tion algorithm exists for BBQP with dominance ratio larger than 1 — 2* ° -(m+n) ^ un i ess p_]\rp 
We then analyze some powerful local search algorithms and show that they can get trapped at 
\Q • a local maximum with objective function value less than a f. One of our approximation algo- 

rithms has an interesting rounding property which provides a data dependent lower bound on 
the optimal objective function value. A new integer programming formulation of BBQP is also 
£T} ' given and computational results with our rounding algorithms are reported. 

o ' 

Keywords: quadratic programming, boolean variables, heuristics, worst-case analysis, domi- 
nation analysis. 

1 Introduction 



The bipartite boolean quadratic programming problem (BBQP) is to 

Maximize f(x, y) = x T Qy + cx + dy 
subject to x G {0, l} m ,y G {0, l} n 

where Q = (qij) is an m x n matrix, c = (c\, C2, • • • , c m ) is a row vector in R m , and d 
{d\ , d2, ■ ■ ■ , d n ) is a row vector in R n . Without loss of generality, we assume that m < n. 
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BBQP has applications in data mining, clustering and bioinformatics |35j . approximating a 
matrix by a rank-one binary matrix |10l I32j . mining discrete patterns in binary data [211 [32], 
solving fundamental graph theoretic optimization problems such as maximum weight biclique [4"[l34j. 
maximum weight cut problem on a bipartite graph [26], maximum weight induced subgraph of a 
bipartite graph [26] , and computing approximations to the cut-norm of a matrix [3] . 

BBQP is closely related to the well-studied boolean quadratic programming problem (BQP) [7[ 
EH [39]: 

Maximize f(x) = x T Q'x + c'x 
subject to x G {0, l} n , 

where Q' is an n x n matrix and c' is a row vector in R n . BBQP can be formulated as a BQP with 
n + m variables [26] and hence the resulting cost matrix will have dimension {n + m) x (n + m). This 
increase in problem size is not desirable especially for large scale problems. On the other hand, we 
can formulate BQP as a BBQP by choosing 

1 1 

Q = Q' + 2MI, c = -c' - Me and d = -c - Me, (1) 

where / is an n x n identity matrix, e £ R n is an all one row vector and M is a very large 
number [26j. Thus, BBQP is a proper generalization of BQP which makes the study of BBQP 
further interesting. An instance of BBQP is completely defined by the matrix Q and vectors c 
and d and hence it is represented by &(Q,c,d). Thus, <$^(Q,0, 0) represents a BBQP with no 
terms cx or dy in the objective function. Such an instance is referred to as a homogeneous BBQP. 
Relationships between BBQP and its homogeneous version are considered in |26| . 

BBQP is trivial if the entries of Q,c and d are either all positive or all negative. BBQP is 
known to be NP-hard |25j since the maximum weight biclique problem (MWBP) is a special case 
of it. Approximation hardness results for MWBP are established by Ambuhl et al. [1] and Tan |34j . 
Performance ratio for approximation algorithms for some special cases of BBQP are discussed by 
Alon and Naor [3] and Raghavendra and Steurer [27]. Results extensive experimental analysis of 
algorithms for BBQP are reported by Karapetyan and Punnen |19] and Glover et al. [13]. Punnen, 
Sripratak, and Karapetyan |26j studied BBQP and identified various polynomially solvable special 
cases. Various classes of valid inequalities and facet defining inequalities for the polytope associated 
with BBQP are obtained by Sripratak and Punnen [33J. 

Worst case analysis of approximation algorithms (heuristics) are carried out normally through 
the measure of performance ratio |37j . Other important measures include differential ratio [8], 
dominance ratio |114 117], dominance number [401 [Tlj. comparison to average value of solutions [5j 
|36"1 [3"Tj [30l [28] etc. Our focus in this paper is on domination analysis and average value based 
analysis of approximation algorithms for BBQP. Berend et al. [7J eloquently argues the importance 
of domination analysis in the study of approximation algorithms. 

Let & be the family of all solutions of BBQP and it is easy to see that |J^| = 2 m+n . The 
average objective function value g/(Q,c,d) of all the solutions of BBQP is given by g/(Q,c,d) = 
2~(m+n) ^ \^gf f(x, y). The idea of comparing a heuristic solution to the average objective func- 
tion value of all the solutions as a measure of heuristic quality for combinatorial optimization 
problems was originated in the Russian literature in the early 1970s. Most of these studies are 
focussed on the traveling salesman problem and the assignment problem (e.g. Rublineckii [28], 
Minina and Perekrest [22] . Vizing [38] . Sarvanov and Doroshko |29[ 130]). In the western literature, 
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Gutin and Yeo [16], Grover |14| . Punnen et al. [25] . Punnen and Kabadi |24j . Deneko and Woeg- 
inger [9] studied the traveling salesman problem and identified heuristics that guarantee a solution 
with objective function value no worse than the average value of all tours. Such a solution has 
interesting domination properties and hence the approach is also relevant in dominance analysis 
of heuristics. For recent developments on domination analysis, we refer to the excellent research 
papers [TJ CEH [15]. Gutin and Yeo [16], Sarvanov [51] . and Angel et al. [5] studied heuristics 
for the quadratic assignment problem with performance guarantee in terms of average value of 
solutions. Similar analysis for the three-dimensional assignment problem was considered by Sar- 
vanov [31], for the Maximum clique problem by Bendall and Margot [6], and for the satisfiability 
problem by Twitto [36]. Berend et al. [7] considered dominance analysis by including infeasible 
solutions. Other problems studied from the point of view of dominance analysis and average value 
based analysis include graph bipartition, variations of maximum clique and independent set prob- 
lems (5j E] and the subset-sum problem [7] . For information on dominance results and linkages 
with the development of heuristic algorithms based on very large scale neighborhood search, we 
refer to [T]. 

A solution with objective function value no worse than the average value of a solution with 
high probability can be obtained by repeated random sampling. However, it should be pointed out 
that even algorithms that performs well in practice could produce solutions with objective function 
value inferior to the average value of a solution. We observed this particularly in the case of the 
BBQP. Thus, a worst case performance of a heuristic that guarantees a solution with objective 
function value no worse than the average value of a solution is a useful measure to be included 
when studying worst case behavior of heuristic algorithms for combinatorial optimization problems. 

Let (x, y), (x°, y°) € JF. Then (x°,y°) dominates (x,y) if f(x,y) < f(x°,y°). Let T be a 
heuristic algorithm for BBQP that produces a solution (x T ,y r ). Define W r = {(x,y) G & : 
f{x, y) < /(x r , y r )}. Let / be the collection of all instances of BBQP. Then the dominance number 
and dominance ratio of V are defined respectively as 



The concept of dominance ratio in the analysis of heuristics was proposed by Glover and Pun- 
nen [llj. Prior to this work, Zemel [40J considered different measures to analyze heuristic algo- 
rithms one of which is equivalent to the dominance number. Hassin and Kuller |17j also considered 
similar measures in analyzing heuristic algorithms. 

In this paper we obtain a closed form formula to compute s/(Q,c,d) in 0(mn) time. We also 
show that any solution to BBQP with objective function value no less than sf (Q, c, d) dominates 
2 m+n ~ 2 solutions. Such a solution is called no worse than average solution. Two algorithms of 
complexity 0(mn) are developed to compute no- worse than average solutions. Thus, the dominance 
ratio of these algorithms is at least One of these algorithms have interesting rounding property 
which provides data dependent lower bounds. The problem of computing a solution with objective 
function value no worse than the median of the objective function values of all solutions is shown 
to be NP-hard. Further, we show that, unless P=NP, for any fixed rational number a > 1, 
no polynomial time approximation algorithm exists for BBQP with dominance ratio larger than 

(l-q) , . \ 

1 - 2 ^ m+n >. We also analyze some very powerful local search algorithms and show that, in worst 
case, such algorithms could get trapped at a locally optimal solution with objective function value 
less than s/ (Q, c, d). Finally we provide a new integer programming formulation of BBQP and the 
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resulting LP relaxation solution could be used to initiate our rounding algorithms. Computational 
results are also provided using the rounding algorithms which establish that the algorithms are 
good candidates to obtain very fast starting solutions in complex metaheuristic algorithms. 

Through out this paper, we use the following notations and naming conventions. We denote 
= {1, 2, . . . , m} and J/ = {1, 2, . . . , n}. The ith component of a vector is represented simply 
by adding the subscript i to the name of the vector. For example, the ith component of the vector 
x* is x*. The set {0, l} n is denoted by B n and [0, l] n is denoted by U n for any positive integer n. 
For any positive integer m, an m- vector of all l's is denoted by l m and an m- vector of all 0's is 
denoted by m . 



2 Average value of solutions and dominance properties 

Note that there are 2 m candidate solutions for x and 2™ candidate solutions for y. Then the 
solutions in the family & can be enumerated as & = {(x k ,y ) : k = 1, 2, . . . , 2 m , £ = 1,2, . . . , 2 n }. 
The next theorem gives a closed form expression to compute £/(Q,c,d) in 0(mn) time. 

Theorem 1. s/(Q, c, d) = - ^ E *j + \ E Ci + \ E d i ' 
Proof. Let rj = 2 m and v = 2 n . Then 

k=l 1=1 

= ^EE EE + E + E d ^ 

k=i e.=i \i^Jt jr'e./K je./K y 

= E E^E^E^ + E c *EE^ + E^EE^ 

\iG^f JG./K fe=l £=1 £=1 k=l je,A" k=l 1=1 



2m+n 



^E E^ + ^E^ + 2 1 E^- 



□ 



If either x = or y = 0, then x T Qy = and such an (x, y) is called a trivial solution. Note that 
f[x,y) need not be equal to zero for trivial solutions. All remaining solutions are called nontrivial 
solutions. Maximizing fix, y) over trivial solutions is straightforward and thus one can restrict 
attention to non-trivial solutions only. The number of nontrivial solutions is (2 m — l)(2 n — 1). Let 
A(Q,c,d) denote the average value of all nontrivial solutions for £?(Q,c,d). 
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i2"^i / 2 ra 1 2 m 1 \ 

Corollary 2. A(Q, c, d) = ———— £ £ Qij + -_ £ c , + __ £ ^ . F U r- 

2"i+« 

tfcer, A(Q,0"\O n ) = -— - —— — — ^(Q,0 m ,0«) and lim lim A(Q,c,d) = #/(Q,c,d). 
(2 — 1)(2 — 1) n->oom-4oo 

From Corollary [2j the asymptotic behavior of A(Q,c,d) and g/(Q,c,d) are the same. Thus, 
hereafter we focus our attention on s/(Q, c, d) only. 

Let <S = {(x,y) : x e {0, l} m ,y € {0, l} n , /(x, y) < *f(Q,c,d)}. Thus, S? consists of all 
solutions of BBQP that are no worse than average. 

Theorem 3. |Sf| > 2 m+n - 2 . 

Proof. For any solution (x, y) G & , let P(x, y) = {(x, y), (x, l n — y), (l m — x, y), (l m — x, l n — y)}. 
It can be verified that P(x, y) = P(x, 1™ — y) = P(l m — x, y) = P(l m — x,l n — y) and P(x, y) / 
P(x',y') if x' $l {x, l m — x} or y' ^ {y, l n — y}. Thus, we can partition the solution space & into 
i 2 m+n = 2 m+n-2 dis j oint sets p^fc^fc^ J~ = 1,2, . . . , 2 m+n - 2 = oj, say. Note that 

f(x, y) + f(l m -x,y) + f{x, l n -y) + f(l m - x, l n - y) = £ £ ?y + 2 £ c, + 2 £ ^ 

JG^jG^ iG^ jG.yC 

= 4^(Q,c,d). (2) 

From equation ([2]), it follows immediately that 

min {f(x, y), /(l m - x, y),f(x, l n - y),f(l m - x,l n - y)} < *f(Q, c, d). 

Thus, from each P{x k ,y k ), k = 1, 2, ... u, choose a solution with smallest objective function value 
to form the set D\. By construction, f(x,y) < g/(Q,c,d) for all (x,y) € D\. Since \D\\ = 2 m+n ~ 2 , 
the result follows. □ 

The lower bound on <S established in Theorem [3] is tight. To see this, consider the matrix Q 
defined by 
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— 1 if i = m, j = n, 
otherwise. 



and choose c and d as zero vectors in R m and R n , respectively. Then £/(Q,c,d) = ^ and the 
set of solutions (x,y) with f(x,y) < g/(Q,c,d) is precisely = {(x,y)\x m = y n = 1}. Clearly, 
|£f | = 2 m+n ~ 2 and hence the bound obtained in Theorem [3] is the best possible. 

Theorem [3] establishes that any algorithm that guarantees a solution with objective function 
value no worse than £?(Q,c,d) dominates 2 m+n " 2 solutions of the BBQP &>(Q,c,d). In other 
words, the domination ratio of such an algorithm is at least 1/4. Our next theorem establishes 
an upper bound on the dominance number of any polynomial time approximation algorithm for 
BBQP. 

Theorem 4. Unless P=NP, no polynomial time algorithm for BBQP can have dominance number 

. I m-j-n I 

more than 2 m+n — 2 L a J for any fixed rational number a > 1. 
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Proof. We show that a polynomial time algorithm f2 for BBQP with dominance number at least 
2 m+n — 2 L ^ J + 1 can be used to compute an optimal solution to BBQP. Without loss of generality, 
assume a = | where a and b are positive relatively prime integers with a > b. Consider an instance 
^(Q, c, d) of BBQP. Let Q* = (q*-) be an abm x abn matrix where 



if i £ M and j € JV ', 
otherwise. 



Likewise, let c* and d* be vectors in R abm and R such that 

I otherwise I otherwise. 

It is easy to verify that from any optimal solution to the BBQP instance ^(Q*, c* , d*) an optimal 
solution to &(Q, c, d) can be recovered. The total number of solutions of 0*(Q*,c*,d*) is 2 abi - m+n ^ 
of which at least 2 ab ( m+n ) - ( m+n ) are optimal. So the maximum number of non-optimal solutions 

ig 2 ab(m+n) _ 2 ab(m+n)-(m+n) Solve thg gggp ingtance <^(Q* )C * )C f) Using ft and let (s*,jf) 

be the resulting solution. By hypothesis, the objective function value of (x*,y*) is not worse 

ab(m+n) „ 

than that of at least 2 ab ( m+n ) - 2 »/» + 1 = 2 ab ( m+n ) - 2 fe ( m+n ) + 1 solutions. Since a > 6, 
we have 2 afe ( m+n ) - 2 b2 ( m+n ) + 1 > 2 ab ( m+n ) - 2 ab ( m + n M m + n ) . Thus, (x*,y*) must be optimal 
for &(Q* ,c* ,d*). From (x*,y*), an optimal solution to &(Q,c,d) can be recovered by simply 
taking the first m components of x* and first n components of y*. The result now follows from 
NP-completeness of BBQP. □ 

Theorem H] implies that unless P=NP, no polynomial time approximation algorithm for BBQP 

(l-q) / , \ 

can have dominance ratio more than 1 — 2 <* ( m + n ) f or ari y nxec l rational number a > 1. 

Although we have a closed form formula for computing the average value of all solutions to 
BBQP, we now show that computing the median value of all solutions is NP-hard. 

Since |J^"| is even, there are two values of median, say 9\ and #2, where 9\ < 62- A median 
finding algorithm could simply produce Q\ or 62 , but we may not know precisely, the output is 
either 9\ or 62 ■ 

Theorem 5. Computing a median of the objective function values of BBQP is NP-hard. 

Proof. Suppose we have a polynomial time algorithm to compute a median of the objective function 
values of BBQP. We will show that this algorithm can be used to solve the PARTITION problem, 
which is defined as follows: Given n positive integers a%, 02, . . . , a n , determine if there exists a 
partition Si and 52 of iV = {1,2,..., re} such that YljeSi a i = Sjgs 2 a .r From an instance of 
PARTITION, construct an instance of BBQP as follows: Choose c as the zero vector, dj = aj 
for j = 1,2,..., re. Choose ^# = {1,2}. Define qij = aje and q2j = — Oj-e, where e is a very 
small positive number. For each subset H of N, let y H be its characteristic vector, i.e. y H € B ri 
and y? = 1 if and only if j G H. For each choice of H, we can associate four choices for x as 
x = (0, 0), x = (1, 0), x = (0, 1) or x = (1, 1). Thus, for each H, we get the following solutions Fh = 
{((0,0), y*), ((1,0), y H ), ((0,1), y H ), ((1,1), y H )}. Now, /((0,0),y*) = E jeH <h>f(0-, 0),O = 

(1 + e) E j&H «i> /((0, 1), V H ) = (1 - e) Ej- 6 h «i an d /((!> 1). J/ H ) = E, e n a r Th ^s, /((0, 1), y H ) < 
/((0,0),y H ) = /((l,l),y H ) < f((l,0),y H ). There are 2™ choices for H and hence there are 2 n+2 
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different solutions for the BBQP constructed. For each subset H of N, let g(H) = YljeH a i an< ^ 
G = {g{H) : H C N}. We first observe that G has two medians and these median values are the 
same and equal to \ dj if and only if N has the required partition. This follows from the 

fact that for any H C N, either g(H) < \ J2 je N a j ^ d( N \ H ) or 9( N \ H ) <l J2 je N a j ^ d( H )- 
Let u\ < o"2 < • • • < cr K be an ascending arrangement of distinct g(H), H C N and let 
W k = {H C iV : cj(iT) = Ofc}. Note that |W fc | = |W K+ i_fc| and a k + a K+ i_ fc = Y.j&N a j- Thus, 
the required partition exists if and only if median of {a\, a^, . . . , <r K } = 5 Sjew a i- Consider an 
ascending arrangement of f(x,y) for all solutions (x,y) of the BBQP constructed. This can be 
grouped as blocks of values B\ < B% < ■ ■ ■ < B K where the block B k has the structure 

repeated \Wk\ times repeated \ W^\ times 

, * v , * V 

cjfe(l-e) = crft(l-e) = • • • = crfc(l-e) < a k = a k = ■ ■ ■ = <Jk < cr k (l+e) = ^(l+e) = • • • = a k (l+e) 

s v y 

repeated 2 1 Vl/fe | times 

for k = 1,2, ... ,k. Thus, median of {a±, 01, . . . , a K } is the same as median of the objective function 
values of BBQP. Thus, the required partition exists if and only both the median values are the same 
and equal to | X^jgTV a r The proof now follows from the NP-completeness of PARTITION. □ 

It may be noted that the above theorem does not rule out the possibility of a polynomial time 
algorithm with dominance number 2 m+n_1 . 



3 Average value of solutions and local search 

In this section we consider two natural local search heuristics for BBQP and show that the solution 
produced could have objective function value worse than si (Q, c, d). One of the popular heuristics 
for BBQP is the alternating algorithm proposed by many authors [21 ( 110 1 H9]. The algorithm starts 
with a candidate solution x° and try to choose an optimal y°. Then fix y° and tries to find the 
best candidate for x, say x 1 yielding a solution (x , y°). These operations can be carried out using 
the formulas 

1 if HjVj + c » > 0, 
otherwise, 

Now fix x = x 1 and try to choose the best y = y 1 and the process is continued until no 
improvement is possible by fixing either x variables or y variables. The algorithm terminates 
when a locally optimal solution is reached. From experimental analysis [19], it is known that this 
algorithm produces reasonably good solutions on average. To the best of our knowledge, worst-case 
behavior of this algorithm has not been investigated. 

Theorem 6. The objective function value of a locally optimal solution produced by the alternating 
algorithm could be arbitrarily bad and could be worse than £/(Q,c,d). 

Proof. Choose m = n, c = d = 0, set qu = l,q nn = M and qij = for all other combinations of i 
and j. Choose the starting solution x\ = 1 and x® = for i 7^ 1. The algorithm will choose = 1 
and yj=0 for j ^ 1. Now (x°,y°) will be locally optimal solution with objective function value 1, 



1 if ^ qijX? + dj > 0, 



ifj = \ ie-# and x\ 

otherwise, 
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but the optimal objective function value is M + 1 for any M > 0. The average cost of a solution is 
M f 1 . Thus, for M > 3 the solution produced by the alternating algorithm is less than g/(Q,c,d). 
In fact, by choosing M large, the solution can be made arbitrarily bad. □ 

Despite this example, it is easy to see that a solution produced by the alternating algorithm 
dominates at least 2 m + 2 n — 1 solutions. 

Let us now consider a more general neighborhood, which is a variation of the ^-exchange 
neighborhood studied for various combinatorial optimization problems. For any (x°,y°) G J^", let 
N hk be the set of solutions in & obtained by flipping at most h components of x° and at most 
k components of y°. If h = m or k = n, we ignore "," in the definition of N h,k . Note that 

|N^| = (Ej=o(?)) (E!U(?))» and the best solution in this neighborhood can be identified 
in polynomial time for fixed h and k. A more powerful neighborhood is N a = N mQ U N an and 

m = 2 m E J =o Q + 2 n E « =q (m) _ e « =q ^ e « =q ( n) Again) thig neighborhood can algo be 

searched for an improving solution in polynomial time for fixed a [19] . It may be noted that a 
solution produced by the alternating algorithm is locally optimal with respect to the neighborhood 
N° = N m0 UN 0n . Glover et al. [[3] considered the neighborhoods N 1 , N 2 , and N 1 ' 1 . They provided 
fast and efficient algorithms for exploring these neighborhoods supported by detailed computational 
analysis. They also considered tabu search algorithms using these neighborhoods in a hybrid 
form. Computational results with these algorithms provided very high quality solutions, improving 
several benchmark instances. Nonetheless, our next theorem shows that even such very powerful 
local search algorithms could provide solutions with objective function values that are inferior to 
), c, d) even if we allow a to be a function of n. 



Theorem 7. A locally optimal solution to BBQP with respect to the neighborhood N a = N ma UN' 
could be worse than average for any a < [^J . 

Proof. Consider the matrix Q defined as 



an 



A if i = m,j = n, 

— 1 if i = m or j = n but ^ (m, n), (3) 
a otherwise. 



and choose c and d as zero vectors. Without loss of generality, we assume m = n. Otherwise, 
we can extend the matrix Q into an n x n matrix by adding n — m rows of zeros and extending 
the vector c into an n- vector by making the last n — m entries zeros. We also assume that n is a 
multiple of 5. Consider the solution (x°,y°) where = y^ = 1 and all other components are zero. 
Also, let a = ^ and assume n > 6. 

Let iV£(0) be the set of all x € B n with x n = obtained by flipping exactly r entries of x° and 
iV£(l) be the set of all x € B n with x n = 1 obtained by flipping exactly r entries of x . Define Ny(0) 
and N r y (0) analogously. Note that for any (x, y) G N an , x G iVJ(O) U JVJ(l) and y G JVJ(O) U JV^(l) 
for some < r < a, < s < n and for any (a?, y) G N nQ , x G JVJ(0) U iVJ(l) and y G iV'(O) U JV"(1) 
for some < r < n,0 < s < a. Thus, for (x, y) G iV a we have 
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f(x,y) 



(r-l)(s-l)a if x G N^(0),y G A^ s (0) for (r, s) € 7£ U /" 

(ra-l)(s-l) if x G iV£(l),y G A^ s (0) for (r, s) G J£U/£, 

(r-l)(so-l) if x G A£(0),i/ G A^ s (l) for (r, s) G /£Ui£, 

rsa-r-s + A if x G A£(l),l/ G N°(l) for (r, s) G 1° U I™. 



where Ig 
„o „,o 



{0,1,..., p} x {0,1,..., <?}. 



Thus, (x , y ) is locally optimal with respect to N a if and only if 

(r-l)(s-l)a < A for all (r,a) G /" U l£ 
(ra-l)(s-l) < A for all (r,s) G U 1%, 
(r-l)(so-l) < A for all (r, s) G 1° U l£, 
rsa-r-s< for all (r, s) G 1° U IJJ. 



It is not very difficult to verify that conditions 



(El), ([7]), and (JU), are satisfied if 



(q - l)(n - l)a < A, 
(n - l)(aa - 1) < A, 
(q — l)(no — 1) < A, and 
aan — a — n < 0. 



(4) 



(5) 
(6) 
(7) 
(8) 



(9) 
(10) 

(11) 
(12) 



Choose a = ^. Then inequality (112]) holds and inequality ([9]) implies inequalities ([10]) and (pT 
Choose A = (a — l)(n — 1)^. Then (x°,y°) is locally optimal. Now, 



£/(Q, c, d) - f{x°, y°) = \{{n- l) 2 a — (2n — 2) + A) — A 

1 /2 o 58 

= — (-n 2 H n - 12) > for n > 6. 

4n 5 5 

This completes the proof. □ 

As an immediate corollary, we have the following result. 

Corollary 8. For any fixed h and k, the objective function value of a locally optimal solution with 
respect to the neighborhood N could be worse than £/(Q,c,d) for sufficiently large m and n. 

Proof. Choose a = max{/i, k}. Then a locally optimal solution with respect to N a is not worse 
than a locally optimal with respect to N ftfc . The result now follows from Theorem [71 □ 

These examples motivates us to develop polynomial time algorithms for BBQP that guarantee 
a solution with objective function value no worse than £/ (Q, c, d). 
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4 Algorithms with no worse than average guarantee 



We first consider a very simple algorithm to compute a solution with objective function value 
guaranteed to be no worse than £/(Q,c,d). The algorithm simply takes fractional vectors x G U n 
and y G U n and applies a rounding scheme to produce a solution for BBQP. Let x G U m and 
y G V n . Extending the definition of f(x,y) for 0-1 vectors, define 

f(x,y) = YY lij^yj + Y CiXi + Y d 3 y r 

Consider the solutions y* G B™ and x* G B m given by 

1 if dj + V] qijXi > 0, 



Vn 



otherwise, 



(13) 



1 if Ci + ^2 qijy* > 0, 
otherwise. 



(14) 



Note that x* is the optimal 0-1 vector when y is fixed at y* , and equation (I13D rounds the y to y* 
using a prescribed rounding criterion. The process of constructing (x*,y*) from (x,y) thus called 
a round-y optimize-x algorithm or RyOx-algorithm. The next theorem establishes a lower bound 
on the objective function value of the solution produced by the RyOx-algorithm. 

Theorem 9. f(x*,y*) > f{x,y). 

Proof. 



fix,y) = YY 1v x iVj + Y CiXi + Y d i y i 



i£.4 



Vi + Y 



Y ( Y qi i Xi + d i ) 

-SIS qi i Xi + d i ) y *i + Y ° 



C^Xi 



(by construction of y* 



Y d ^ Y\Y <w*i + 



Ci I X?' 



^ + Y [Y + c * 

= f(x*,y*). 



(by construction of x*) 



□ 



Note that (x*,y*) can be constructed in 0(mn) time whenever x and y are rational numbers. 
We can also round x first to obtain x° G M m and choose optimal y = y° by fixing x at x°. This is 
done using the rounding scheme given by the following equations: 
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T 



1 if Ci + ^ lijVj > 0, ( 1 if dj + ^ > 0, 

je./K and y ° = < ie„# 

otherwise, [0 otherwise, 

Theorem 10. f(x°,y°) > f(x,y). 

The proof of Theorem [TH] follows along the same line as Theorem [9] and hence omitted. The 
process of constructing (x°,y°) is called round-x optimize-y algorithm or RxOy-algorithm. The 
complexity of the RxOy-algorithm is also 0{mn). 

Corollary 11. A solution (x,y) for BBQP satisfying f(x,y) > ^(Q,c,d) can be obtained in 
0(mn) time. 

Proof. Let xi = 1/2 for all i E ^# and yj = 1/2 for all j E JV . Then it can be verified that 
f(x, y) = g/ (Q, c, d), but (x, y) is not feasible for BBQP. Now, choose (x, y) as the output of either 
the RyOx-algorithm or the RxOy-algorithm. The result now follows from Theorems l9l or [TUl □ 

In view of Corollary [IT] and Theorem [9j the dominance ratio of the RyOx-algorithm and the 
RxOy-algorithm is at least j. By choosing an appropriate starting solution, we can establish 
improved dominance ratio for these algorithms. 

We now discuss an unexpected upper bound on g/ (Q, c, d). As a consequence, we have yet an- 
other simple scheme to compute a solution that is not worse than average. Let a = ^ ^2j €y y Qij, 
P = J2 ie ^ c i and 7 = Eje.^ dj- 

Theorem 12. <s/(Q, c, d) < max{a + /3 + 7, /3, 7, 0}. 

Proof. Let u and v be real numbers in [0, 1]. Choose x E B m , y E B n be such that Xi = u for 
all i E ^# and yj = v for all j E JV . Then f(x,y) = auv + /3u + >yv = n(u,v), say. Note that 
77(1/2, 1/2) = g/(Q,c,d). Thus, max{r/(u, v) : (u,v) E U 2 } > £/(Q,c,d). Since n(u,v) is bilinear, 
its maximum is attained at an extreme point of the square U 2 . Since these extreme points are 
precisely (0,0), (1,0), (0,1), (1,1), the result follows . □ 

Corollary 13. One of the solutions (l rn , l n ), (1"\ n ), (0 m , l n ), (0 m , n ) of BBQP have an objec- 
tive function value no worse than &/(Q,c,d). 

The proof of this corollary follows directly from Theorem [T2j We can compute a, f3 and 7 in 
0(mn) time and hence we have a solution no worse than average in 0(mn) time. Interestingly, 
if a, f3 and 7 are given, then we can identify a solution to BBQP with objective function value 
no worse than &/(Q,c,d) in O(l) time. The solution produced by Corollary 1131 is trivial and may 
not be of much practical value. Nevertheless, the simple upper bound on (Q, c, d) established by 
Theorem [12] is very interesting and have interesting consequences as discussed below. 

Recall that the alternating algorithm starts with a solution (x°,y°), fix x° and find the best y, 
say y . Then fix y at y 1 and compute the optimal x and so on. Since we initiate the algorithm 
by fixing x first, we call this the x-first alternating algorithm. We can also start the algorithm by 
fixing y° first and the resulting variation is called the y -first alternating algorithm. 

Theorem 14. The best solution amongst the solutions produced by the x-first alternating algorithm 
with starting solutions (l m , l n ) and (0 m ,0 n ) dominates 2 m+n ~ 2 + 3 (2™" 1 ) solutions. 
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Proof. Let (x*,y*) be the best solution obtained. When the starting solution is (l m , 1"), f(x*,y*) > 
max{/(l m , 1"), f(l m , 0")}. Likewise, when the starting solution is (0 m ,0 n ), we have f(x*,y*) > 
max{/(0 m ,0 ri ),/(0 m ,l n )}. Thus, by Theorem QJ/(a;*, y*) > £?(Q,c,d) and hence by Theorem 
(x*,y*) dominates at least 2 m+n ~ 2 solutions. To account for the remaining solutions that are dom- 
inated by (x*,y*) we proceed as follows. 

Construct the set D\ as in the proof of Theorem [3j Recall that \D\\ = 2 m+n ~ 2 . Let D 2 = 
{(l m ,y),(0 m ,y) : y G B n }. Now, \D 2 \ = 2 n+1 . For any y G B n , we have P(0 m ,y) = P{l m ,y) and 
P(0 m ,y) = {(0 m ,y),(l m ,y),(0 m ,l n -y),(l m ,l n -y)}. Thus, D 2 can be partitioned into sets of the 
form P(0 m ,y k ) and there are j2 n+1 = 2 n ~ 1 such sets. Exactly one element from each P(0 m ,y k ) is 
mD x . Thus, \DiDD 2 \ = 2 n ~ l and hence \DiUD 2 \ = 2 m + n " 2 + 2 n+1 -2™- 1 = 2 n ~ 1 (2 m ~ l + 3). □ 

From Theorem [14] the dominance ratio of the x-flrst alternating algorithm is at least | + 2 „f +1 
when starting solutions are selected carefully and the algorithm is applied twice. Note that to 
achieve this dominance ratio, we simply need to perform only one iteration each, when the algo- 
rithm starts with (l m , l n ) and (0 m , ra ). Thus, we can achieve this dominance ratio in polynomial 
time. A similar result can be derived for y-first alternating algorithm and local search algorithm 
with neighborhood N a for any a > 0. Further, similar dominance ratio can be achieved by RxOy- 
algorithm or RyOx-algorithm applied twice, once starting with (l m ,l n ) and then starting with 
(0 m ,0 n ) and choosing the best solution. 

The problem BBQP can be formulated as integer linear programming problem [26] as follows: 

ILP1: Maximize Qij z ij + °i x i + djUj 

Subject to Zij — Xi < 0, i € j £ Jf 
Zij — yj < 0,i G , j € JV 
Zij - Xi- yj > -1, G S, 

Xi G {0, 1} for i G J( , yj G {0, 1} for j G Jf , Zij G {0, 1} for G x JV, 

where S = {ij : qij < 0}. Let (z',x',y') be an optimal solution to the LP relaxation, where 
z' is an m x n matrix with (i,j)th entry z'-, and x' is an m- vector with ith entry x[ and 
y' is an n- vector with jth entry y'-. Recall that for x' G [0, l] m and y' G [0, l] n , 4>(x' ,y') = 
}l, . '/-/•'•!//' + T,ie^ c i x 'i + Ysjejsdjy'j and h(z',x',y') is the optimal objective function 

value of the linear programming relaxation of ILP1. 

For any x G [0, l] m and y G [0, l] n , the solution (x , y 1 ) obtained by the RxOy (RyOx) algorithm 
satisfies cj)(x,y) < f(x 1 ,y l ) < h(z' , x' ,y'). This follows from Theorem [9] and the property of LP 
relaxations. Alternative integer programming formulations of BBQP could provide different LP 
relaxation solutions and hence the resulting RxOy (RyOx) rounding solutions could be different. 
We give below a new integer programming formulation of BBQP by increasing the number of 
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variables. 

ILP2: Maximize ^ ^ + %■ - \wij - + ^ + ^ 

Subject to Uij + 2tijj = Xi + for i G j G jV 
Uij + Ujj < 1, for i G j G =/K 
— 0y + Wij = Xi — yj for i G j G jV 
Zij + Wij < 1 for i € , j E «/K 
Xj G {0, 1} for i G J{ , yj G {0, 1} for j G =yf 
Uij, Vij, Wy, G {0,1} for G x 

Experimental analysis of RxOy (RyOx) algorithms starting from the LP relaxations of ILP1 
and ILP2 are discussed in the next section. 



5 Computational results 

Although the primary focus of the paper is on theoretical analysis of approximation algorithms, we 
have conducted preliminary experimental analysis with the RxOy and RyOx rounding algorithms 
to examine features that are not clear from theoretical analysis. The alternating algorithm and local 
search algorithms are thoroughly analyzed from an experimental analysis point of view in [T9J and [?] 
and hence are not considered in our experimental study. The rounding algorithms being simple 
constructive heuristics, they are not expected to outperform more sophisticated algorithms that 
employ powerful neighborhoods within a metaheuristic framework. Nonetheless, it is interesting to 
explore the behavior of these algorithms in the light its theoretical properties and the potential for 
computing initial solutions for more advanced solution improvement algorithms. 

The data set used in our experiments are smallsize instances of random problems, biclique 
instances, matrix factorization instances, and maxcut instances from [19]. The algorithms are 
implemented in C# and tested on a DELL PC with Windows 7 operating system, Intel i7 processor 
and 16GB of memory. All CPU times reported are in milliseconds and do not include input-output 
times. We used the solution obtained by the linear programming relaxation of ILP1 as the fractional 
vector (x, y) to initiate the rounding process. 

Our first set of experiments were aimed to identify the percentage of problems where LP re- 
laxation produced optimal solutions. This data is important since for such problems, rounding 
(heuristic algorithms) is irrelevant. The results of these experiments are summarized in Figure 1. 
The trend shows that as the problem size increases, the number of problems where the LP relax- 
ation produced optimal solution decreases. This decrease is more rapid in some class of problems 
compared to random instances. 

Let us now discuss the experimental results with the RxOy and RyOx algorithms. The details 
are summarized in Tables 1 to 5. In each table, the column "best" reports the best known objective 
function value of the problem. Most of these values are optimal. The column "LP obj" contains the 
upper bound obtained by the LP relaxation. The columns "xy" and "yx" respectively contains the 
objective function value of the solution produced by the RxOy algorithm and the RyOx algorithms. 
It may be noted that the solution produced by these algorithms are close to the best known solution 
values and the running time is very minimal. This makes these algorithms good candidates to be 
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50 100 150 200 50 100 150 200 



ran mn 
Figure 1: Random Figure 2: Matrix Factorization 

used to generate starting solutions for more sophisticated algorithms. The table also provide insight 
into the value s$ (Q, c, d) in comparison to the best known objective function value. The column 
"avg" reports srf (Q, c, d) and the column "Avg+" reports the value of the upper bound on stf (Q, c, d) 
provided by Theorem [T2l Note that Avg+ is the objective function value of the solution reported 
in Corollary 1131 These are certainly inferior solutions. The column ll XiUj" represents f(x,y) for 
the LP relaxation solution (x,y). The experimental results conclusively demonstrate the power of 
RxOy and RyOx algorithms as a very fast stand alone heuristics or as an algorithm for generating 
starting solutions in more complex algorithms. 



Objective Time, ms 



Instance 


Best 


LP obj. 


x iUj 




xy 


Avg+ 


Avg 


LP 




xy 


20 x 50 


13555 


21179 


437.5 


12715 


12208 


2160 


624.0 


151 


171 


146 


25 x 50 


13207 


25405 


759.0 


12226 


11016 


2911 


505.0 


69 


77 


87 


30 x 50 


15854 


29713 


634.0 


13992 


14824 


3219 


778.0 


102 


105 


102 


35 x 50 


14136 


33471 


-2146.0 


12554 


12099 


835 


-2028.5 


137 


143 


131 


40 x 50 


18778 


39444 


42.8 


16562 


15573 


624 


125.8 


168 


183 


170 


45 x 50 


22057 


44760 


-401.0 


17688 


19660 


610 


-440.0 


212 


236 


200 


50 x 50 


23801 


50390 


1576.0 


21645 


22178 


5115 


1576.0 


234 


272 


233 



Average 17341 34909 128.9 15340 15365 2211 162.9 154 169 153 



Table 1: Random 
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Objective 



Time, ms 



Instance 


Best 


LP obj. 


b & J 


y x 


xy 


Avg+ 


Avg 


LP 


y x 


x y 


20 x 50 


18341 


29999 


-1911363.0 


4933 








-1911363.0 


47 


49 


47 


25 x 50 


24937 


38513 


-3265058.5 


2764 








-3265058.5 


66 


90 


65 


30 x 50 


27887 


51726 


-2981174.8 


9050 


15138 





-2981174.8 


97 


95 


95 


35 x 50 


32515 


57302 


-3782013.0 


3270 


10502 





-3782013.0 


135 


134 


134 


40 x 50 


33027 


61966 


-4562391.5 











-4562391.5 


133 


131 


133 


45 x 50 


37774 


67923 


-5789225.5 


4420 








-5789225.5 


140 


141 


140 


50 x 50 


30124 


78745 


-5004969.5 


5517 


5833 





-5004969.5 


200 


195 


200 


Average 


29229 


55168 


-3899456.5 


4279 


4496 





-3899456.5 


117 


119 


1160 



Table 2: Max Biclique 



Objective Time, ms 



Instance 


Best 


LP obj. 


XiVj 


yx 


xy 


Avg+ 


Avg 


LP 


yx 


xy 


20 x 50 


6983 


10704 


-315.3 


5223 


5481 





-393.3 


110 


109 


111 


25 x 50 


8275 


13866 


-678.5 


7630 


7350 





-678.5 


169 


166 


164 


30 x 50 


10227 


18958 


-140.0 


8378 


8685 





-140.0 


232 


234 


228 


35 x 50 


11897 


20590 


-543.8 


11071 


9008 





-543.8 


283 


283 


277 


40 x 50 


14459 


23110 


892.3 


13216 


13433 


3569 


892.3 


506 


524 


498 


45 x 50 


13247 


24393 


-984.5 


12496 


9513 





-984.5 


430 


431 


428 


50 x 50 


15900 


28875 


356.3 


14075 


15357 


1425 


356.3 


628 


622 


621 


Average 


11570 


20071 


-201.9 


10298 


9832 


713 


-213.1 


337 


338 


333 



Table 3: Max Induced Subgraph 



Objective Time, ms 



Instance 


Best 


LP obj. 




yx 


xy 


Avg+ 


Avg 


LP 


yx 


xy 


20 x 50 


9008 


21266 


-791.0 


3838 


4928 





-780.0 


93 


83 


89 


25 x 50 


10180 


27546 


-1205.0 


4062 


6410 





-1352.0 


123 


130 


120 


30 x 50 


13592 


37676 


-1172.0 


5808 


7560 





-274.0 


205 


184 


183 


35 x 50 


14024 


40900 


-1196.0 


5816 


8578 





-1105.0 


331 


267 


269 


40 x 50 


17610 


45948 


1133.0 


9424 


12082 


3568 


1784.0 


334 


290 


297 


45 x 50 


15252 


48492 


-1209.0 


6204 


9790 





-1956.0 


441 


402 


402 


50 x 50 


19580 


57408 


723.0 


12332 


10364 


1446 


723.0 


591 


538 


536 


Average 


14178 


39891 


-531.0 


6783 


8530 


716 


-422.9 


303 


271 


271 



Table 4: MaxCut 
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Objective 



Time, ms 



Instance 


Best 


LP obj. 






xy 


Avg+ 


Avg 


LP 


yx 


xy 


20 x 50 


114 


249 


-1.5 


106 


94 





-1.5 


37 


32 


31 


25 x 50 


127 


314 


1.5 


109 


96 


6 


1.5 


45 


40 


42 


30 x 50 


148 


373 


-2.5 


124 


140 





-2.5 


65 


61 


61 


35 x 50 


139 


425 


-13.0 


106 


122 





-13.0 


80 


77 


78 


40 x 50 


210 


512 


12.0 


182 


194 


48 


12.0 


148 


73 


71 


45 x 50 


191 


559 


-3.5 


165 


179 





-3.5 


100 


93 


94 


50 x 50 


217 


637 


11.5 


181 


191 


46 


11.5 


120 


119 


119 


Average 


164 


438 


0.6 


139 


145 


14 


0.6 


85 


71 


71 



Table 5: Matrix Factorization 
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We have also conducted experiments with RxOy and RyOx rounding algorithms using the LP 
relaxation solution of ILP2. For random instances, the LP relaxation solution of ILP1 and ILP2 
were different but for all other test instances, they produced same solutions most of the time. 
In some random instances considered, the rounding algorithms produced better solutions when 
started ILP2 LP relaxation solution in comparison to the solutions produced from ILP1 LP relax- 
ation. However, in general we did not see significant performance difference and hence preferred 
ILP1 formulation, considering its smaller size. 

Note that the RxOy and RyOx rounding algorithms can be initiated using any x° £ U m and 
yO ^ jjjn an j no £ necessar iiy a fractional solution for an LP relaxation. For i = 1,2, ... ,m let 
7i = Ci + Ylje^rQij an d f° r 3 = 1,2, ... ,n let Sj = dj + ^2 ie ^qij- Also, ran[a,b] represents a 
uniformly distributed random number in the interval [a, b]. Using these values we generate random 
vectors (x°,y°) as follows: 
Type 1 random vectors: Choose 

_ Jran(0, 0.5] if 7i < 0, _ Jran(0, 0.5] if Sj < 0, 

— \ and y j — \ 

lran(0.5, 1] otherwise, I ran(0.5, 1] otherwise. 

Type 2 random vectors: These are generated using weighted random numbers with weight 
proportional to the corresponding 7 j or Sj values. Let p\ = max{| 7 j| : i £ M} and p2 = max{|<5j| : 
j e JV}. Now choose 



^0.5 - ^ran(0,0.5] if 7i < 0, 
7i Pi and 
— ran(0.5, 1] otherwise, 

L Pi 



0.5 - ^-ran(0, 0.5] if Sj < 0, 
5- 

-j-ran(0.5, 1] otherwise. 

P2 



The RxOy and RyOx algorithms using the above choices of starting solution (x°,y°) can also 
be used to construct starting solutions for advanced algorithms. The built-in randomness generates 
good solutions that can be embedded in metaheristics with multiple starts. Systematic experimental 
analysis of such sophisticated algorithms is beyond the scope of this paper. 



6 Conclusion 

In this paper we studied approximation algorithms for BBQP which is a generalization of the well 
known BQP. Various approximation algorithms are analyzed using averaged value based measures 
and domination analysis. It is demonstrated that very powerful local search algorithms could get 
trapped at poor quality local minimum even if we allow exponential time in searching a very large 
scale neighborhood. Some of the proof techniques used arc simple yet elegant and could be of 
use in domination analysis of heuristics for other related problems. Experimental results with two 
construction algorithms are also given. A natural question for further investigation is to close the 
gap between non-approximability bounds and lower bounds on domination ratio. Since BBQP is 
not as well studied as BQP, there are many other avenues for further investigation and we are 
currently investigating further properties of BBQP. 
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